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ABSTRACT: The metaphor of holey adaptive land- 
scapes provides a pictorial representation of the process of 
speciation as a consequence of genetic divergence. In this 
metaphor, biological populations diverge along connected 
clusters of well-fit genotypes in a multidimensional adap- 
tive landscape and become reproductively isolated species 
when they come to be on opposite sides of a "hole" in the 
adaptive landscape. No crossing of any adaptive valleys is 
required. I formulate and study a series of simple models 
describing the dynamics of speciation on holey adaptive 
landscapes driven by mutation and random genetic drift. 
Unlike most previous models that concentrate only on 
some stages of speciation, the models studied here de- 
scribe the complete process of speciation from initiation 
until completion. The evolutionary factors included are 
selection (reproductive isolation), random genetic drift, 
mutation, recombination, and migration. In these mod- 
els, pre- and post-mating reproductive isolation is a conse- 
quence of cumulative genetic change. I study possibilities 
for speciation according to allopatric, parapatric, peri- 
patric and vicariance scenarios. The analytic theory sat- 
isfactorily matches results of individual-based simulations 
reported by Gavrilets et al. (1998). It is demonstrated 
that rapid speciation including simultaneous emergence 
of several new species is a plausible outcome of the evo- 
lutionary dynamics of subdivided populations. I consider 
effects of population size, population subdivision, and lo- 
cal adaptation on the dynamics of speciation. I briefly 
discuss some implications of the dynamics on holey adap- 
tive landscapes for molecular evolution. 

Speciation has traditionally been considered to be one 
of the most important and intriguing processes of evolu- 
tion. In spite of this consensus and significant advances 
in both experimental and theoretical studies of evolution, 
understanding speciation still remains a major challenge 
(Mayr 1982a; Coyne 1992). The main reason for such 



a discouraging situation is that direct experimental ap- 
proaches, which are widely used for solving other prob- 
lems of evolutionary biology, are not effective for studying 
speciation because of the time scale involved. Experi- 
mental work necessarily concentrates on distinct parts of 
the process of speciation intensifying and simplifying the 
factors under study (Rice and Hostert 1993; Templeton 
1996). In situations where direct experimental studies are 
difficult or impossible, mathematical modeling has proved 
to be indispensable for providing a unifying framework. 
Although numerous attempts to model parts of the pro- 
cess of speciation have been made, a quantitative theory 
of the dynamics of speciation is still missing. Currently, 
verbal theories of speciation are far more advanced than 
mathematical foundations. As often is the case with ver- 
bal theories (both scientific: and otherwise), different de- 
duced (or induced) aspects of speciation are emphasized 
by different workers resulting in confusion and contro- 
versy. The situation is not helped by the absence of gen- 
eral agreement on a species definition (e.g., Claridge et 
al. 1997). 

Here, I attempt to develop some foundations for a gen- 
eral dynamical theory of speciation. One possible ap- 
proach to this goal would be to begin with a species defini- 
tion, then to define speciation accordingly and to develop 
an appropriate dynamical model. I do not think such ap- 
proach would be very useful, due to a lack of generality. 
My models are not based on a specific "species concept". 
I reason that species are "different" with respect to some 
characteristics, and that whatever these differences, they 
have a genetic basis. Thus, modeling the dynamics of spe- 
ciation is equivalent to modeling the dynamics of genetic 
divergence. I use a bottom-up approach: begin with a 
model incorporating a range of factors thought to lead to 
speciation (e.g. selection, mutation, population subdivi- 
sion etc.) and then try to interpret its dynamic behavior 
in terms of different "species concepts". As expected, 
many aspects of speciation that are emphasized by dif- 
ferent species concepts (such as reproductive isolation, 
separate genotypic clusters or common evolutionary tra- 
jectories) emerge from the same processes. This clearly 
indicates that different species concepts are not mutually 
exclusive. 

The choice of a modeling approach depends upon the 
purpose of the model. A common view in (evolution- 
ary) biology is that mathematical models are mainly use- 
ful for making predictions that can be used in experi- 
mental work. Although such a pragmatic approach is 
probably what should be expected in contemporary so- 
ciety, a model's testable predictions are not necessarily 
its main contribution to science. Insights provided by 
models, their ability to train our intuition about complex 
phenomena, to provide a framework for studying such 
phenomena and to identify key components in complex 
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systems are at least as important as specific predictions. 
For these purposes the most useful tools are simple mod- 
els and metaphors. 

Sewall Wright's (1932) metaphor of "rugged adaptive 
landscapes" is a well-known and widely used metaphor 
in evolutionary biology. In the standard interpretation, 
a rugged adaptive landscape is a surface in a multidi- 
mensional space that represents the mean fitness of the 
population as a function of gamete (or allele) frequencies 
which characterize the population state (see Fig. la). It is 
envisioned that this surface has many peaks and valleys 
corresponding to different adaptive and maladaptive pop- 
ulation states, respectively. The population is imagined 
as a point on the surface which is driven by selection up 
hill but can get stuck on a local peak. Two general points 
about scientific metaphors should be kept in mind. The 
first is that specific metaphors (as well as mathematical 
models) are good for specific purposes only. The second is 
that accepting a specific metaphor necessarily influences 
and defines the questions that are considered to be im- 
portant. The metaphor of "rugged adaptive landscapes" 
is very useful for thinking about adaptation. However, 
its utility for understanding speciation is questionable. 
From a pragmatic point of view, the process of splitting a 
population into two different species is impossible to de- 
scribe in a framework where a population is the smallest 
unit. Finer resolution is necessary for describing the split- 
ting of populations. Accepting the metaphor of rugged 
adaptive landscapes immediately leads to a problem to 
be solved: how can a population evolve from one adap- 
tive peak to another across an adaptive valley when selec- 
tion opposes any changes away from the current adaptive 
peak? Wright's solution to this problem, his shifting- 
balance theory (Wright, 1931, 1982), does not seem to 
be satisfactory (Gavrilets 1996; Coyne et al. 1997). 
Provinc (1986), Barton and Rouhani (1987), Whitlock 
et al. (1995), Gavrilets (1997a) and Coyne et al. (1997) 
discuss other weaknesses of Wright's metaphor. I have 
argued elsewhere (Gavrilets 1997a) that his metaphor of 
rugged adaptive landscapes with its emphasis on adap- 
tive peaks and valleys is, to a large degree, a reflection of 
the three-dimensional world we live in. Both genotypes 
and phcnotypcs of biological organisms differ in numerous 
characteristics, and, thus, the dimension of "real" adap- 
tive landscapes is much larger than three. Properties of 
multidimensional adaptive landscapes are very different 
from those of low dimension. Consequently, it may be 
misleading to use three-dimensional analogies implicit in 
the metaphor of rugged adaptive landscapes in a multi- 
dimensional context. I believe that understanding speci- 
ation requires a different metaphor. 

The metaphor of "holey" adaptive landscapes. An in- 
dividual organism can be considered as a combination of 
genes. All possible combinations of genes form a geno- 
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Figure 1: Adaptive landscapes, a. A rugged 
adaptive landscape. The main emphasis is on 
"peaks" and "valleys" corresponding to different 
well-adapted and maladaptive genotypes, b. A 
holey adaptive landscape. The main emphasis 
is on clusters of well-fit genotypes that extend 
throughout the genotype space. All other geno- 
types are treated as "holes". 

type space (which, mathematically, can be represented 
by a hypercube). In discussing the evolution of popula- 
tions, it is useful to visualize each individual as a point in 
this genotype space. Accordingly, a population will be a 
cloud of points, and different populations (or species) will 
be represented by different clouds. Selection, mutation, 
recombination, random drift and other factors change the 
size, location and structure of these clouds. To construct 
an adaptive landscape one assigns "fitness" to each geno- 
type (or each pair of genotypes) in genotype space. Dif- 
ferent forms of selection and reproductive isolation can 
be treated within this conceptual framework. For exam- 
ple, fitness can be a genotype's viability (in the case of 
viability selection); it can be fertility or the probability 
of successful mating between a pair of genotypes (in the 
case of fertility selection or premating isolation, respec- 
tively). A finite population subject to mutation is likely 
to be represented by genotypes with fitnesses within a 
fitness band determined by the balance of mutation, se- 
lection and random drift. A general property of adap- 
tive landscapes with a very large number of dimensions is 
that genotypes with fitnesses within a specified band form 
connected "clusters" that extend throughout the geno- 
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type space (Gavrilets 1997; Gavrilets and Gravner 1997). 
Thus, populations can evolve and diverge along bands of 
highly-fit genotypes without going across the states with 
a large number of low- fit genotypes (that is without cross- 
ing any adaptive valleys). 

The metaphor of "holey" adaptive landscapes puts spe- 
cial emphasis on these clusters of highly- fit genotypes, dis- 
regarding fitness differences between them and treating 
all other genotypes as "holes" (Gavrilets 1997; Gavrilets 
and Gravner 1997). The justification for the latter is a 
belief that selection will be effective in moving the popu- 
lation away from these areas of genotype space on a time 
scale that is much faster than the time scale for speciation. 
Accordingly, microevolution and local adaptation can be 
viewed as the climbing of the population from a "hole" to- 
wards the holey adaptive landscape, whereas macroevolu- 
tion can be viewed as a movement of the population along 
the holey landscape with speciation taking place when the 
diverging populations come to be on opposite sides of a 
"hole" in the adaptive landscape. In this scenario, there is 
no need to cross any "adaptive valleys" ; reproductive iso- 
lation between populations evolves as an inevitable side 
effect of accumulating different mutations. The metaphor 
of holey adaptive landscapes can be traced to a verbal 
two-locus two-allele model of reproductive isolation pro- 
posed by Dobzhansky (1937) and similar ideas discussed 
by Bateson (cited in Orr, 1997), Muller (1942), Maynard 
Smith (1970, 1983), Nei (1976), Barton and Charlesworth 
(1984), Kondrashov and Mina (1986). For more discus- 
sion of this metaphor see Gavrilets (1997ab) and Gavrilets 
and Gravner (1997). Orr (1995) and Gavrilets (1997a) 
gives a summary of relevant experimental evidence. The 
metaphor of "holey" adaptive landscapes is illustrated 
graphically in Figure lb. 

Mathematical models for holey adaptive landscapes. 
Here I briefly review previously published work on the 
evolutionary dynamics on "holey" adaptive landscapes. 
Nei (1976) and Wills (1977) were the first to present for- 
mal analyses of the Dobzhansky model. Nei et al. (1983) 
studied one- and two-locus multi-allele models with step- 
wise mutations and considered both postmating and prc- 
mating reproductive isolation. In their models geno- 
types were reproductively isolated if they were different 
by more than 1 or 2 mutational steps. In these situations, 
speciation was very slow. They conjectured, however, 
that increasing the number of loci may significantly in- 
crease the rate of speciation. Bengtsson and Christiansen 
(1983) presented a deterministic analysis of mutation- 
selection balance in the Dobzhansky model. Bengts- 
son (1985), Barton and Bengtsson (1986) and Gavrilets 
(1997b) analyzed the properties of hybrid zones arising 
under Dobzhansky-type epistatic selection. Wagner et al. 
(1994) considered a two-locus, two-allele model of stabi- 
lizing selection acting on an epistatic character. For a 



specific set of parameters, the interaction of epistasis in 
the trait and stabilizing selection on the trait resulted 
in a fitness "ridge". The existence of this ridge sim- 
plified stochastic transitions between alternative equilib- 
ria. Gavrilets and Hastings (1996) formulated a series 
of two- and three-locus Dobzhansky-type viability selec- 
tion models as well as models for selection on polygenic 
characters. They studied these models in the context of 
founder effect speciation and noticed that the existence 
of ridges in the adaptive landscape made stochastic di- 
vergence much more plausible. Similar conclusions were 
reached by Gavrilets and Boake (1998) who studied the 
effects of premating reproductive isolation on the plau- 
sibility of founder effect speciation. Higgs and Derrida 
(1991, 1992) proposed a model where the probability of 
mating between two haploid individuals is a decreasing 
function of the proportion of loci at which they are differ- 
ent. Here, any two sufficiently different genotypes can be 
considered as sitting on opposite sides of a hole in a ho- 
ley adaptive landscape. These authors as well as Manzo 
and Peliti (1994) studied this model numerically assum- 
ing that the number of the loci is infinite, the loci are 
unlinked and highly mutable, and mating is preferential. 
Orr (1995) and Orr and Orr (1996) studied speciation 
in a series of models in which viability of a diploid or- 
ganism depends on the number of heterozygous loci. All 
these papers postulated the existence of ridges of highly- 
fit genotypes. Gavrilets and Gravner (1997) studied a 
general class of multilocus selection models and showed 
the existence of ridges to be inevitable under fairly gen- 
eral conditions. Independently, a similar conclusion was 
reached in Reidys et al. (1997). Most previous studies of 
the dynamics of speciation on holey adaptive landscapes 
were numerical. To develop a dynamical theory of speci- 
ation it is desirable to have a simple model that can be 
treated analytically. 

The Model 

I consider finite populations of haploid individuals with 
discrete, non-overlapping generations. I assume that re- 
production involves gene exchange (amphimixis) between 
individuals. The restriction to haploids is for algebraic 
simplicity. Models for diploids will be discussed later. In- 
dividuals are different with respect to L possibly linked 
diallelic loci. Without any loss of generality each individ- 
ual's genotype can be represented as a sequence of 0's and 
l's. Let l a = (I", ...,12) where is equal to or 1, be such 
a sequence for an individual a. In standard population 
genetics models, the population state is usually described 
in terms of gamete frequencies. In systems with many 
loci such an approach is not practical. For instance, with 
10 diallelic loci there arc 2 10 different gametes. Thus, one 
would need to analyze more than 1000 coupled equations. 
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Another complication follows from the fact that even in 
very large populations with hundreds of thousands indi- 
viduals each specific genotype is represented only by a 
small number of copies or is not represented at all. Thus, 
the notion of a gamete frequency in multilocus evolution 
might be very difficult to justify. Here, I will be inter- 
ested in the levels of genetic variation within subpopu- 
lations and genetic divergence between subpopulations. 
Both can be characterized in terms of genetic distance 
d defined as the number of loci at which two individu- 
als are different. More formally, the genetic distance d a $ 
between individuals a and j3 is 

L 

d a/3 = J2( l i- l i) 2 - (!) 

8=1 

Genetic distance d is the standard Hamming distance. 
It is analogous to the number of segregating sites in a 
sample of two gametes, which is widely used in molecu- 
lar evolutionary genetics (Li, 1997), and to the number 
of heterozygous loci in a diploid organism. Genetic dis- 
tance d is also closely related to the notion of the overlap 
q between two sequences, d = ^(1 — q), which is com- 
monly used in statistical physics (e.g., Derrida and Peliti 
1991). I model the expected dynamics of average genetic 
distances within and between populations, using D w for 
the former and Df, for the latter. 

I assume that reproductive isolation is caused by cumu- 
lative genetic change. I will use a very simple symmetric 
model which is closely related to the models discussed 
above and which allows one to treat both pre- and post- 
mating isolation within the same framework. I posit that 
an encounter of two individuals can result in viable and 
fecund offspring only if the individuals are different at 
no more than K loci. Otherwise the individuals do not 
mate (premating reproductive isolation) or these offspring 
are inviable or sterile (postmating reproductive isolation) . 
More formally, I assign "fitness" w to each pair of individ- 
uals depending on the genetic distance d between them 

f „ / 1 for d < K, . . 

w ( d ) = { terd>K. ® 

(see Appendix for an outline of more complicated ap- 
proaches). In this formulation, any two genotypes differ- 
ent at more than K loci can be conceptualized as sitting 
on opposite sides of a hole in a holey adaptive landscape 
(cf. Higgs and Derrida 1991, 1992). At the same time, a 
population can evolve to any reproductively isolated state 
by a chain of single locus substitutions. The adaptive 
landscape corresponding to this model is both "holey" 
and "correlated" . The latter means that the probability 
that two genotypes are reproductively isolated correlates 
with the genetic distance between them. In Nei et al. 



(1983) and Gavrilets and Boake (1998) models individ- 
uals separated by more than one mutational step were 
reproductively isolated which corresponds to K = 1. The 
neutral case (no reproductive isolation) corresponds to K 
equal to the number of loci. 

The mathematical model presented above was inter- 
preted as describing sexual haploid populations with fit- 
nesses assigned to pairs of individuals depending on the 
genetic distance between them. However, there is an al- 
ternative interpretation in that the model describes ran- 
domly mating diploid populations. In the diploid case, 
the genetic distance @) between the two gametes forming 
an individual is equivalent to individual's heterozygosity, 
and fitness function (Q) specifies fitness as a function of 
individual heterozygosity. Therefore, most conclusions of 
this paper will also be applicable to situations when post- 
mating reproductive isolation is in the form of reduced (or 
zero) viability of hybrids due to incompatibility of the 
genes they receive from their parents (Wu an Palopoli 
1994). 

Dynamics in the neutral case 

Before developing a theory for the dynamics of specia- 
tion in the above model, it is illuminating to start with 
the neutral case. Here I summarize some relevant results 
that are presented in (or follow directly from) classical 
papers (Watterson 1975; Li 1976; Slatkin 1987; Strobeck 
1987). Let fi be the probability of mutation per locus per 
generation. The approximations below assume that mat- 
ing is random, the number of loci L is large, but (i is very 
small so that the probability of mutation per individual 
per generation v = L\i << 1. The migration rate m and 
the inverse of the population size 1 /N are small as well. 

Genetic variation within an isolated population. Let 
us consider an isolated population of size Nt- The ex- 
pected change in the average genetic distance within the 
population per generation is 

AD w =2v-^, (3) 

where the first term in the right-hand side is the contribu- 
tion of mutation whereas the second term is the random 
drift reduction of D w . Asymptotically, a mutation-drift 
equilibrium is reached with 

D* w = 9 = 2vN T . (4) 

Genetic divergence between isolated populations. Let 
us consider several isolated populations of arbitrary size. 
The probability that a specific mutation gets fixed in a 
population is one over the population size. Different mu- 
tations will get fixed in different populations resulting in 
their genetic divergence. The average genetic distance 
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between any two of them increases with the rate equal 
twice the mutation rate per gamete 



AD,, 



2v. 



(5) 



The rate of neutral divergence does not depend on pop- 
ulation sizes. In particular, it is the same independent of 
whether there are many small populations or a few large 
populations. Because the number of loci L is finite, an 
indefinite increase of D w , which is implied by equation 
(^|) is impossible. This equation as well as equation ([|) 
above and equations (||) below approximate the dynamics 
when genetic distances D w and D b are small relative to 
the number of loci L. To treat the general case, one has to 
substitute v for v(l — 2D W /L) in equations (0) and ( |6a|) 
and for v(l — 2D b /L) in equations (||) and (|6b|). With 
a finite number of loci, the genetic distance D b between 
isolated populations approaches L/2 asymptotically. 

Subdivided populations. The effect of migration on the 
average genetic distances depends on the spatial structure 
of populations. Assume that a population of size Nt is 
subdivided into n subpopulations of size N = Nt/ti and 
that a proportion m > of individuals migrate to any of 
the other n — 1 subpopulations. The expected changes 
in the average genetic distances within and between sub- 
populations are 



AD W = 2v + 2m(D b - D u 
2m 



EL 
N ' 



AD b = 2v 



n- 1 



(D w - D b ) 



(6a) 
(6b) 



Equations (|^) assume that v,m and l/N are small. 
Asymptotically, a mutation-migration-drift equilibrium is 
reached with 



D* 

n 



i) 



2v 



(7a) 
(7b) 



where 9 is given by equation (^). The average genetic 
distance within a subpopulation (of size N) does not de- 
pend on the number of subpopulations n and migration 
rate m and is the same as is expected in a single popula- 
tion with size Nt- The average genetic distance between 
sub-populations increases with the population subdivi- 
sion and decreasing migration. Figure 2 illustrates the 
dynamics of neutral divergence in a system of two sub- 
populations. 

Peripheral population. Assume a "peripheral" popu- 
lation of size N is receiving migrants from a very large 
"main" population. Genetic variation in the main popu- 
lation is assumed to be constant (and not influenced by 
migration from the peripheral population) . The expected 
changes in the average genetic distances within the pe- 




Figure 2: Dynamics of D w and Db in the neutral case. 
Population size N = 100 (solid lines) and N = 200 (dashed 
lines). The rate of migration is m = 0.01, the mutation rate 
per individual is v = .0384. Initially, D w — Db = 0. 



riphcral population, D Wl and between the peripheral and 
main populations, D b are 

AD w = 2v + 2m(D b -D w )-^, (8a) 



AD b = v + m(D Q - D b 



(8b) 



where Do is the average genetic distance within the main 
population and m is the proportion of individuals in the 
peripheral population replaced by migrants from the main 
population. Asymptotically, a mutation-migration-drift 
equilibrium is reached with 



„ _ 2Nm 
Dw = l + 2Nm D(h 

D* b =D + -, 



(9a) 
(9b) 



where the former equation assumes that genetic variation 
in the main population is sufficiently large (Do >> vN) 
and the number of migrants, Nm, is not too small. The 
average genetic distance within the peripheral population 
is always larger than that for an isolated population of 
its size (D^ > 2vN) . If the number of migrants is large 
(Nm >> 1), the average genetic distance within the pe- 
ripheral population is about the same as in the "main" 
population. 

Dynamics with reproductive isolation 

The main feature of the model for reproductive isola- 
tion introduced above and other models of holey adaptive 
landscapes is the existence of chains of equally-fit combi- 
nations of genes separated by single substitutions that 
extend throughout the genotype space. These chains can 
be though of as "neutral paths" in the adaptive land- 
scape. It is important to realize, however, that the exis- 
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tence of "holes" in a holey adaptive landscape makes the 
actual dynamics of genetic divergence not neutral. In this 
section I summarize some analytical results on the evo- 
lutionary dynamics in the case of reproductive isolation 
described by equation (|^). Gavrilets et al. (1998) have 
studied the possibilities for speciation in this model nu- 
merically. Details of the analytical methods used are out- 
lined in the Appendix. To derive the dynamic equations 
below, I have used the same assumptions as described 
above at the beginning of the section on the neural case 
substituting the assumption of random mating for the 
assumption of random encounters. In addition, I have 
assumed that the distributions of genetic distances both 
within and between populations are Poisson. There are 
several sets of approximations resulting in Poisson distri- 
bution of genetic distances. In the present contest, the 
weakest set seems to be the assumption that genetic vari- 
ation at each locus is small most of the time (rare-alleles 
approximation) and that the population is approximately 
at linkage equilibrium. These assumptions are standard 
in analyzing the dynamics of multilocus systems under 
the joint action of selection, mutation, and random drift 
(e.g., Barton 1986; Barton and Turclli 1987; Burger et al. 
1989; Gavrilets and de Jong 1993). The fit of individual- 
based simulations with analytic predictions is satisfactory 
both at qualitative and quantitative levels (see below and 
Gavrilets et al. 1998). 

Genetic variation within an isolated population. After 
the population becomes polymorphic at K loci, new mu- 
tations are selected against when rare because individuals 
carrying them have a reduced probability of producing vi- 
able and fecund offspring. Selection experienced by indi- 
vidual loci underlying reproductive isolation is frequency- 
dependent (and is similar to that arising in the case of 
underdominant selection on a diploid locus). The change 
in D w per generation in an isolated population of size N 
is approximately 



where 



AD W = —sD„, 



s = 



2v- 



N ' 



D 



K 



T(K + 1,D W ) 



(10) 



(11) 



and T(x,y) is an incomplete gamma function (e.g., 
Gradshteyn and Ryzhik, 1994). The value of D* w at 
the mutation-drift-selection equilibrium can be found by 
equating the right-hand side of ( |l0| ) with zero and solving 
for D w . Figure 3a illustrates the dependence of D w on 
the parameters of the model. This figure indicates that 
the equilibrium values of D w are close to the correspond- 
ing neutral predictions (Q) if K is larger than 2 or 3 times 
(where = 2Nv is the average genetic distance within a 
finite population in the neutral case). Figure 3b gives the 



values of the effective selection coefficient s. With moder- 
ately large K (that is with K > 10), s is very small. The 
effective selection coefficient s can also be thought of as 
the strength of induced selection on each locus underlying 
reproductive isolation. Figure 3b shows that very strong 
selection on the whole genotype (implied by the existence 
of complete reproductive isolation at finite values of K) 
results in very weak selection at the level of individual 
loci. 





Figure 3: (a) Average genetic distance D w maintained by 
mutation-selection-drift balance in an isolated population of 
size N as a function of K for v = 0.0384. The circles, squares, 
diamonds and triangles give estimates from individual-based 
simulations for N = 100, 200, 400 and 800, respectively (thirty 
runs for each parameter configuration), (b) Effective selection 
coefficient s in the case of infinite population size for two values 
of v. (c) The genetic load 1 — w for parameters values as in 
Figure b. 

The mean fitness of the population, W Wl can be defined 
as the proportion of pairs of individuals that can mate 
and produce fertile and viable offspring (cf., Nei et al. 
1983). For a population with an average genetic distance 

1,D W ) 



w w — 



T(K+1) 



(12) 



where T(x + 1) is a gamma function (e.g., Gradshteyn 
and Ryzhik, 1994. For integer x, T(x + 1) = xl). Figure 
3c shows that in spite of relatively high levels of genetic 
variation maintained in the population, the genetic load 
(that is the proportion of reproductively isolated pairs 
of individuals, 1 — w w ) is very low. This seems to be a 
general property of holey adaptive landscapes (cf. Wills 
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1977; Bengtsson and Christiansen 1983). 

Genetic divergence between isolated populations. Even 
after the genetic distance within an isolated population 
has reached an equilibrium level, the population keeps 
evolving as different mutations get fixed. As a conse- 
quence, isolated populations will continuously diverge ge- 
netically. The asymptotic rate of divergence of two iso- 
lated populations of size N each is 



AD b = 2vR, 



where 



R 



2e~ s VS 
H erf{sfS) 



(13a) 



(13b) 



is the rate of divergence relative to the neutral case. Here, 
S = Ns/2, s is defined by equation ( pi] ) with D w corre- 
sponding to the mutation-selection-drift equilibrium, and 
erf(x) is the error function (— 2/^/tt J Q exp(—y 2 )dy). In 
the neutral case, s = 0, U = 1/N and equation ( |13a| ) re- 
duces to equation (||) . Figure 4 illustrates the dependence 
of the relative rate of divergence R on model parame- 
ters. In the neutral case, the rate of genetic divergence 
ADb does not depend on the population size (equation 
^|). In contrast, with reproductive isolation the rate of 
divergence decreases with increasing population size. Af- 
ter the population becomes polymorphic at K loci, new 
mutations are selected against when rare. Genetic drift 
operating in finite populations overcomes the effect of se- 
lection and allows genetic divergence. For example with 
K = 20 and v = .0384, a population of size N = 800 will 
accumulate about 5 substitutions per 1000 generations. 
A few thousand generations will be sufficient for Db to 
exceed K significantly. In contrast, very large randomly 
mating populations will diverge very slowly. 

Figure 4 indicates that the rate of substitutions is close 
to the corresponding neutral predictions if K is larger 
than 2-3 times 9 (9 — 2Nv). Note that as in the neutral 
case considered above, an implicit assumption in equation 
( |j"3| ) is that genetic distance D b is small relative to the 
number of loci L. In the general case, AD b = 2v(l — 
2D b /L)R and D b approaches L/2 asymptotically. 

At what moment can the two diverging population be 
considered as two different species? The answer obviously 
depends on what one means by a species. Let us say that 
the two populations are different species if the proportion, 
Wb, of encounters between individuals from different pop- 
ulations that can result in mating and viable and fertile 
offspring is less than a small number 7. (This definition 
uses the biological species concept.) During initial stages 
of divergence, this proportion can be approximated by 
the right-hand side of equation ([l2]) with Db taking the 
place of D w . Figure 5 shows the minimum genetic dis- 
tance between populations required for speciation as a 
function of K for several values of 7. One can see that a 




Figure 4: The rate of difergence R relative to the neu- 
tral case in an isolated population of size TV as a function of 
K for v — .0384. The circles, squares, diamonds and tri- 
angles give estimates from individual-based simulations for 
N = 100, 200, 400 and 800, respectively (thirty runs for each 
parameter configuration) . 

genetic distance between the populations on the order of 
2 or 3 times K will be sufficient for the status of separate 
"biological" species. Note that there is very little effect 
of the magnitude of 7. 

Speciation in a subdivided population. In the determin- 
istic limit (that is with Nt —> 00), the genetic variation 
of a subdivided population can be maintained by migra- 
tion. This can happen if initially alternative alleles are 
close to fixation in different subpopulations and selection 
is sufficiently strong relative to migration (e.g., Karlin 
and McGregor 1972). Let k be the number of loci at 
which alternative alleles are close to fixation in different 
subpopulations. Respectively, L — k will be the number 
of loci at which the same allele is close to fixation in dif- 
ferent subpopulations. In the deterministic limit, k does 
not change. In the n-island model the dynamics of D w 
and Db are described by equations 



AD W = -sD w + 2v + 2m e (D b - D u 

2m e 



(14a) 



AD b = -s(D b - k) + 2v + UD W - D b ), (14b) 

n — 1 

where s is defined by equation (pi]), and the "effective" 
migration rate 

m e = m =r - (15) 

w w 

if k < K and m e — otherwise. Here, w w is given by 
equation (12) above whereas w b = T(K + 1 — k,Db — 
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Figure 5: Minimum genetic distance between populations 
for speciation for 7 = .0001 , .001 and .01 (solid lines from top 
to bottom). Also shown is the diagonal D — K (dashed line) 



k) /T{K + 1 — k) is the probability that two randomly cho- 
sen individuals from different populations are not repro- 
ductively isolated. The effective migration rate m e can 
be thought of as half the probability of mating between 
individuals from different subpopulations. With no re- 
productive isolation (with very large K) or no genetic di- 
vergence between subpopulations (with Db w D w , k = 0), 
the effective migration rate is equal to the actual migra- 
tion rate (m e — m). Comparing equations ( p^ ) with their 
neutral analogs (|J) shows that reproductive isolation re- 
sults in two effects. First, it directly reduces genetic vari- 
ation within subpopulations and genetic divergence be- 
tween subpopulations. These effects are described by the 
first terms in the right-hand side of equations (|l4j). Also, 
reproductive isolation reduces the gene flow between pop- 
ulations. Given that Db > D w , then m e < m reflecting 
the fact that genes brought by migrants have a smaller 
probability of being incorporated in the resident popula- 
tion. In the deterministic limit, both D w and Db always 
evolve to finite equilibrium values. 

Random genetic drift results in two effects. First, it re- 
duces the genetic variation within sub-populations by the 
amount D w /N. The dynamic equation for D w becomes 

AD W = -sD w + 2v + 2m e (D b -D w )- (16) 

Second, genetic drift might change k. The expected 
change in k can be approximated as 



where R is given by equation ( |13b| ). The first term in 
the right-hand side of equation (|17|) can be thought of as 
the rate at which an allele that is initially rare in both 
sub-populations becomes close to fixation in one of the 
subpopulations. This rate was found by Lande (1979) 
using a diffusion approximation and assuming that mi- 
gration is weak. The second term is the rate at which 
the loci with different alleles initially close to fixation in 
different subpopulations become fixed for the same allele 
in both of them. To find this term I used Barton and 
Rouhani's (1987) method. 

Depending on parameter values and initial conditions 
there are two different dynamic regimes. In the first 
regime, both D w and Df, evolve to finite values, which 
are smaller than those in the neutral case (and which are 
much smaller than the number of loci L). Here, selec- 
tion (reproductive isolation) reduces genetic divergence 
both within and between subpopulations. In the sec- 
ond regime, D w stays small (relative to L) whereas Db 
increases effectively indefinitely (to values order L/2). 
Here, selection (reproductive isolation) reduces the ef- 
fective migration rate to zero resulting in speciation. 
These dynamics can be understood in the following way. 
Changes in D w and Db induced by selection are expected 
to happen on a faster time scale than changes in k in- 
duced by random genetic drift. Thus, D w and Db values 
should be close to the equilibrium values predicted by 



equations (14b, 16) when k is treated as a constant. The 
dynamic behavior depends on whether k reaches a finite 
equilibrium value or keeps increasing. In the latter case, 
the effective migration rate m e reduces to zero, the rate 
of change of k approaches 2vR with Db increasing at the 



Afc = 2vR2 



-Nm. 



2km e R(e/2) 



Nm c 



(17) 



same rate (cf. equation 13b). 

Figure 6 illustrates the dynamics observed by numer- 
ically iterating the model equations. The iterations 
started with all N individuals identical. During the first 
1000 generations there were no restrictions on migra- 
tion between subpopulations and the whole population 
evolved as a single randomly mating unit (cf. Gavrilets 
et al. 1998). The average genetic distance within the 
population D w evolves according to equation (|l^). Start- 
ing with generation 1000, restrictions on migration were 
introduced and the dynamics are described by equations 
(14b- 1?]) afterwards. After generation 1000, each of these 
figures has two sets of three curves corresponding to two 
different values of the parameter(s) under consideration. 
The curves within each set represent Df,, D w and k. With 
migration rate m = 0.01, all these variables evolve to- 
wards finite equilibrium values (see Fig. 6a) whereas with 
a smaller migration rate (m = 0.001), Db and k increase 
effectively indefinitely signifying that speciation has taken 
place. Thus, reducing migration makes speciation more 
plausible. Figure 6b shows that increasing mutation rate 
(from v — 0.0384 to 5 times this value) has a similar ef- 
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feet. These two figures describe the dynamics expected 
in a system of two subpopulations. Figure 6c compares 
the dynamics observed in a population subdivided into 2 
and 4 subpopulations. This figure shows that increasing 
population subdivision makes speciation more plausible. 
Note that the process of genetic divergence described in 
Fig. 6c results in a simultaneous emergence of 4 species. 
In the cases where speciation takes place (as signified by 
continuous increase in the genetic distance between sub- 
populations), the curves representing D b and k are par- 
allel meaning that asymptotically the genetic divergence 
is due to fixation of different mutations in different sub- 
populations. In the cases where speciation does not take 
place, k is close to zero. 

Altogether, at the qualitative level the results presented 
in Fig. 6 correspond to both biological intuition and the 
results of individual-based simulations in Gavrilets et al. 
(1998). At the quantitative level, there is a very good 
fit between simulations and analytical predictions for lev- 
els of genetic variation maintained in subpopulations and 
the asymptotic rate of divergence between subpopula- 
tions. However, the conditions for speciation as predicted 



by iterating equations (14b-17) appear to be more strict 
than those observed in the individual-based simulations 
performed by Gavrilets et al. (1998). For example, for 
parameter values used in Figure 6c, no speciation in a 
system of four subpopulations occurs if m > 0.0035. In 
contrast, in individual-based simulations speciation was 
observed for m = 0.01 (Figure 3b in Gavrilets et al. 1998). 
The main reason for this disperancy seems to be an inad- 



equacy of equation (17) at moderate levels of migration 



(e.g. Lande 1979; Barton and Rouhani 1987). 

Speciation in a peripheral population. Here I consider 
the case of a peripheral population of size N receiving 
migrants from a very large main population. The dynam- 
ics of the average genetic distance within the peripheral 
population, D w , the average genetic distance between the 
peripheral and main populations, D b , and the number of 
diverged loci k are approximated by equations 



AD W = - sD w + 2v + 2m e (D b - D w ) - (18a) 
AD b = -^{D b -k) + v + m e (D - D b ), (18b) 



Afc =vR2' Nm ' - km e R(e/2) 



Nl. 



(18c) 



Here Dq is the average genetic distance within the main 
population. Figure 7 illustrates the dynamics observed 
by numerically iterating equations (|l8|). For Dq I used 
mutation-selection balance values for a very large isolated 
population predicted by equations (|l(], [ll]). The initial 
values of D w and D b were equal to Dq. The parameter 



values in Fig. 7a and Fig. 7b are the same as those that 
resulted in speciation in Fig. 6a and 6b, respectively. The 
outcome of the dynamics is the same - speciation - but the 
rate of divergence is smaller than when all subpopulations 
are uniformly small. This is apparent from the level of 
genetic distance between subpopulations achieved after 
1000 generations of divergence which are about twice as 
small in Figures 7a and 7b relative to those in Figures 6a 
and 6b. 

Discussion 

The theory developed above together with earlier numer- 
ical simulations (see references above) show that rapid 
speciation is a plausible outcome of the evolutionary dy- 
namics in subdivided populations. Here speciation is a 
consequence of two fundamental factors. The first fac- 
tor is the existence of various and possibly significantly 
different highly-fit combinations of genes underlying di- 
verse solutions (genetical, ecological, behavioral, devel- 
opmental etc.) to the problem of survival and reproduc- 
tion. In multidimensional genotype space these combina- 
tions of genes tend to form connected clusters that extend 
throughout genotype space. At the same time these geno- 
types are not mutually compatible - they are separated 
by "holes". The second factor is mutation pressure. Be- 
cause the population size is finite and the number of loci 
is very large whereas the probability of a specific mu- 
tation is very small, different mutations tend to appear 
(and increase in frequency) in different subpopulations 
(cf., Barton 1989; Mani and Clarke 1990). Metaphori- 
cally speaking, mutation tends to tear apart the cloud 
of points representing the population in genotype space. 
Combining genes from two different organisms in one off- 
spring can counteract the disruptive effect of mutation, 
keeping a single randomly mating population together in 
genotype space. But restricting gene exchange as a conse- 
quence of limited migration between subpopulations gives 
mutation a significant advantage. Eventually the popula- 
tion cloud will be broken and smaller clouds representing 
the subpopulations will drift apart in genotype space - 
an event representing speciation. Given sufficient genetic 
divergence, restoring migration to high levels will not re- 
turn the system back to the state of free gene exchange 
between subpopulations which now can be considered as 
different species. It is not necessary to invoke strong selec- 
tion for local adaptation to explain speciation in a subdi- 
vided population, as studied here, or after a founder event 
(Gavrilets and Hastings 1996; Gavrilets and Boake 1998). 
Mutation is ubiquitous. Population size is never infinite 
and, thus, genetic drift is always present. Speciation as 
caused by mutation and random drift should represent 
a null model against which speciation as caused by local 
adaptation can be tested (cf., Nei 1976; Lande 1976). 
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Unlike most previous models that concentrate only on 
some stages of speciation, the model studied here de- 
scribes the complete process of speciation from initiation 
until completion. I assumed that reproductive isolation 
is caused by cumulative genetic change. The model is 
described in terms of dynamic equations for the vari- 
ables analogous to those used in molecular evolution- 
ary biology - the average genetic distances between and 
within subpopulations. Average genetic distances within 
(sub)populations always evolve towards finite equilibrium 
values. Depending on parameter values and initial con- 
ditions average genetic distances between subpopulations 
either converge to a finite equilibrium or increase effec- 
tively indefinitely. The former regime is interpreted as no 
speciation. In the latter regime, three effects take place 
simultaneously: (1) genetic distances between subpopula- 
tions significantly exceed genetic distances within them, 

(2) encounters between individuals from different sub- 
populations do not result in viable and fertile offspring, 

(3) evolutionary changes in a subpopulation do not affect 
other subpopulations. Thus, subpopulations form sepa- 
rate genotypic clusters in genotype space, become repro- 
ductively isolated and undertake changes as evolutionary 
independent units. This regime is interpreted as specia- 
tion according to any of the species concept common in 
the literature (e.g., Mallet 1995; Claridge et al. 1997). 

The dynamic equations derived above describe the ex- 
pected changes in the average genetic distances neglect- 
ing stochastic fluctuations around the expected values. 
The predicted dynamics have two clearly distinct regimes: 
convergence towards a finite equilibrium (no speciation) 
or effectivly indefinite divergence (speciation). Stochas- 
tic fluctuations around the expected values, which are 
present in natural populations (and individual-based sim- 
ulations), make the boundary between these two regimes 
less strict and may result in the population escaping the 
first regime and entering the second regime after some 
time (see Gavrilets et al. 1998). My analysis has been 
based on approximations which are standard in study- 
ing multilocus systems. I assumed that alleles are rare, 
that linkage disequilibrium, mutation and migration rates 
are small and used a theory developed by Lande (1979), 
Walsh (1982) and Barton and Rouhani (1987) for de- 
scribing stochastic transitions driven by random genetic 
drift. The analytic theory presented here fits satisfacto- 
rily with the results of individual-based simulations. The 
model can be used to evaluate qualitative effects of dif- 
ferent factors on the dynamics of speciation, the order of 
magnitude of parameters resulting in or preventing spe- 
ciation, and the time scale involved. According to both 
biological intuition and previous numerical simulations, 
increasing mutation rate and decreasing migration pro- 
mote speciation. Increasing the number of loci has sig- 
nificantly increased the plausibility of speciation relative 



to that in earlier models (Nei et al. 1983; Wagner et al. 
1995; Gavrilets and Hastings 1996). Note that the actual 
number of loci influences the dynamics only through the 
mutation rate per gamete, v, and parameter K . For re- 
alistic parameter values the time scale for speciation can 
be as short as a few thousands or even hundreds of gener- 
ations. This is compatible with rates observed in several 
cases of rapid speciation in natural populations described 
recently (Schluter and McPhail 1992; Yampolsky et al. 
1994; Johnson et al. 1996; McCune 1996, 1997) includ- 
ing the most spectacular case - the origin of hundreds of 
species of Lake Victoria cichlids in 12,000 years (Johnson 
et al. 1996). The model has demonstrated the plausibility 
of speciation with relatively low levels of both initial ge- 
netic variation and new genetic variation introduced into 
the population each generation (both supplied by muta- 
tion). With higher levels of the former (as in laboratory 
experiments on speciation, reviewed by Rice and Hostert 

1993, and Templeton, 1996) or of the latter (for instance 
as a result of natural hybridization, reviewed by Bullini 

1994, Rieseberg 1995, Arnold 1997), the rate of speciation 
is expected to be even higher. 

Local adaptation and speciation 

The model analyzed above shows that rapid speciation in 
a subdivided population can occur even without any dif- 
ferences between selection regimes operating in different 
subpopulations (that is without selection for local adap- 
tation). An important question is how genetic changes 
brought about by selection for local adaptation would af- 
fect the dynamics of speciation (e.g., del Solar 1966; Ayala 
et al. 1974; Kilias et al. 1980; Dodd 1989; Schluter 1996; 
Givnish and Sytsma 1997). These effects will depend on 
whether the genes responsible for local adaptation are 
different from or are the same as the genes underlying 
reproductive isolation. 

Assume first that the two sets of genes are completely 
different. Let the strength of selection per locus induced 
by reproductive isolation be very small so that these loci 
can be considered as effectively neutral. (For the model 
studied here, this seems to be the case if K is larger that 
2-3 times 6, where 9 is the average genetic distance main- 
tained by mutation in a finite population in the neutral 
case.) Then, Birky and Walsh's (1988) results tell us 
that the rate of substitution in these loci will not be af- 
fected by selection on other loci independently of linkage. 
However, given that reproductive isolation is a result of 
genetic incompatibilities, the loci underlying reproduc- 
tive isolation will be under frequency-dependent selection 
against rare alleles which is analogous to underdominant 
selection in diploid populations. Birky and Walsh (1988) 
have shown that linkage to advantageous alleles slightly 
increases the rate of fixation of detrimental mutations. 
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This suggests that selection on linked loci will increase 
the rate of substitutions in the loci underlying reproduc- 
tive isolation and, thus, will promote speciation to some 
degree. No results seem to be known on how linkage to 
advantageous alleles increases the rate of fixation of un- 
dcrdominant mutations or alleles experiencing frequency- 
dependent selection. No quantitative predictions can be 
made here, but most likely if the two sets of loci are not 
extremely tightly linked, effects of selection for local adap- 
tation on the rate of speciation will not be significant. 

Assume now that the loci under consideration 
pleiotropically affect both survival in a given environment 
and reproductive isolation. For instance, this may be 
the case if disruptive selection acts on habitat preferences 
which also define mating patterns (e.g., Rice 1984; Rice 
and Salt 1988) or if the probability of mating between 
individuals depends on the difference in their morpho- 
logical traits which are under direct selection. Let sla 
be the average strength of selection per locus induced by 
selection for local adaptation. Using Walsh's (1982) re- 
sults, the relative rate of fixation of new mutations in an 
isolated population of size N can be approximated (see 
Appendix) as 



R 



7T erf(VS(l + a)) + erf{VS(l- a)) 



(19) 



where S = Ns/2, s is the strength of selection per lo- 
cus induced by reproductive isolation, and a — sla/s. 
With a = 0, equation ( |l9|) reduces to (13b). Figure 8 
illustrates the dependence of R on S and a. Increasing 
a always increases R. Thus, selection for local adapta- 
tion always increases the rate of substitutions and pro- 
motes speciation. With sufficiently strong selection for 
local adaptation (sla > s), the net effect of new alleles 
will be advantageous and their frequencies will tend to 
increase even when rare. In the limit of large popula- 
tion size, the probability of fixation is 2{sla — s). This 
is analogous to the classical results on the probability of 
survival of an advantageous mutant in a very large pop- 
ulation (Haldane 1927; Walsh 1982). The rate of accu- 
mulation of genetic differences will be 2(sla — s)Nv and 
can be significant. Very strong artificial selection for local 
adaptation has been shown to result in rapid evolution of 
reproductive isolation (e.g., del Solar 1966; Kilias et al. 
1980; Dodd 1989). However, the changes brought about 
by moderately strong artificial selection may not exceed 
those resulting from random genetic drift only (e.g. Ringo 
et al. 1985). 

Population subdivision and speciation 

In the models considered here, speciation is a by-product 
of fixation of different alleles in different subpopulations. 



It is well known that the rate of fixation of neutral alleles 
does not depend on population size, that of advantageous 
alleles increases with population size, and that of delete- 
rious or underdominant alleles decreases with population 
size (e.g., Gillespie 1991; Ohta 1992). At the level of 
individual loci, selection induced by reproductive isola- 
tion in the form considered here is similar to underdom- 
inant selection (or frequency-dependent selection against 
rare alleles). Thus, in the absence of selection for local 
adaptation (or with independent loci controlling traits 
for local adaptation) decreasing population size will in- 
crease the rat e of substitutions and promote speciation 
(see equation |l3b| ) . Effects of the population size on the 
plausibility of speciation will be similar even if the same 
loci control both reproductive isolation and locally ben- 
eficial traits given that selection for local adaptation is 
not too strong (sla < s, see equation |l^ and Fig. 9). In 
all these cases, speciation will be driven by mutation and 
random genetic drift and will be fastest if the popula- 
tion is subdivided into small subpopulations. This con- 
clusion about the effect of population subdivision on the 
probability of speciation in Dobzhansky-type models dif- 
fers from that of Orr and Orr (1996). They argued that 
the degree of population subdivision has no effect on the 
rate of speciation if speciation is caused by mutation and 
random drift. Orr and Orr did not consider the actual 
process of fixation of new mutations assuming that it 
will be a simple neutral process. However, the existence 
of holes in the adaptive landscape makes the process of 
substitution non-neutral and new mutations are selected 
against when rare. Such mutations are fixed more easily 
in smaller subpopulations. For the discussion of the ex- 
isting experimental evidence regarding effects of random 
genetic drift on the plausibility of speciation see Rice and 
Hostert (1993) and Templeton (1996). The time scale 
for speciation is short meaning that restrictions on mi- 
gration between subpopulations do not need to be long 
lasting; several hundreds of generations may be sufficient 
for a significant divergence and evolution of reproductive 
isolation. It is quite possible that several new species 
will emerge from a highly subdivided population within 
a short period of time (see above; Gavrilets et al. 1998). 
These theoretical conclusions are consistent with a verbal 
"micro-allopatric" model of speciation suggested for cich- 
lid fishes in the East African Great Lakes (e.g., Reinthal 
and Meyer, 1997). Hoelzer and Melnick (1994) have em- 
phasized that the possibility of simultaneous emergence 
of several new species should be incorporated more ex- 
plicitly in the contemporary methods for reconstructing 
phylogenies. 

If the same loci control both reproductive isolation and 
locally beneficial traits and selection for local adaptation 
is sufficiently strong (sla > s), increasing the population 
size will result in increasing the rate of substitutions (see 
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Fig. 9). In this case, speciation will be driven by selection 
and will be fastest if the population is subdivided into 
a small number (say, two) of large subpopulations (Orr 
and Orr 1996) as implied by the vicariance scenario (e.g. 
Wiley 1988). 

Many species are thought to be represented by a few 
large populations and many smaller "peripheral" popula- 
tions. Mayr (1963, 1982b) proposed the theory of peri- 
patric speciation arguing that speciation is typically ini- 
tiated in small peripheral populations and he attributed 
a special role to genetic drift in this process. Gavrilets 
(1996) has shown that an invasion of a new adaptive com- 
bination of genes is most successful if it is initiated in a 
peripheral population. The results presented here bear 
out Mayr's argument(see Fig. 8). Small peripheral pop- 
ulations will rapidly diverge genetically from the "main" 
large population and speciate. Although differences in se- 
lection regimes between peripheral and main populations 
can accelerate divergence, random genetic drift will be the 
most important factor. On the other hand, if a periph- 
eral population is large enough and is under a selection 
regime that is sufficiently different from the one operating 
in "main" populations, then disruption of gene flow can 
cause evolutionary divergence, perhaps leading to rapid 
speciation, in the absence of contributions from random 
genetic drift (Garcia-Ramos and Kirkpatrick 1997). 

Summarizing, large randomly mating populations will 
diverge genetically and speciate only if there is strong se- 
lection for local adaptation (for instance after a change 
in the environment). In contrast, small populations will 
diverge and speciate even without differences in selec- 
tion regimes between them. Possibilities for speciation 
strongly depend on the geographic structure of the pop- 
ulation. Here, analysis was restricted to the island model 
and the continent-island model. Manzo and Peliti (1994) 
and Gavrilets et al. (1998) present numerical results for 
stepping-stone models. 

Relationship to other speciation models 

Using genetic distance (1) implies the equivalence of loci. 
A general case of non-equivalent loci can be described by 
introducing a [L x L) matrix G = {Gy} of weights and 
defining a generalized distance between individuals a and 

(3 as 

d a P = {l a -l^fGil -I ), (20) 

where l a and l@ are vectors defining the corresponding 
genotypes and superscript T means transpose. Consid- 
ering haploid populations and premating isolation only, 
the model assumes that individuals can mate only if they 
are not too different genetically. Here, the degree of 
reproductive isolation was controlled by cumulative ge- 
netic difference. However, using the generalized distance 
(EOj) allows one to treat models for reproductive isola- 



tion controlled by quantitative traits as well as models 
for sexual selection within the same framework (see Ap- 
pendix). The close relationship between the models of 
speciation as a consequence of "quasi-neutral" divergence 
along ridges in the adaptive landscapes and as a conse- 
quence of sexual selection was already recognized by Bar- 
ton and Charlesworth (1984). 

A fundamental reason for speciation on a holey adap- 
tive landscape is mutation which tends to break the pop- 
ulation into reproductively isolated pieces. Population 
subdivision and the resulting reduction in gene exchange 
facilitates this process. Here, migration rates compati- 
ble with rapid speciation were small (that is speciation 
was allopatric or parapatric). An interesting question is 
whether speciation is possible with much higher migration 
rates. In other words, is sympatric speciation by muta- 
tion and random genetic drift on a holey adaptive land- 
scape possible? Numerical simulations of similar models 
of sympatric speciation where mutation rates were higher 
(Higgs and Derrida 1991, 1992) or the time span studied 
was longer (Wu 1985) or the population size was smaller 
(Gavrilets and Boake 1998) than here, provide an affir- 
mative answer. Adding disruptive selection due to either 
abiotic factors (say, different resources) or biotic factors 
(competition) should create additional pressure on the 
population cloud which might result in rapid sympatric 
speciation. 

Beyond holey landscapes 

Gavrilets and Gravner's (1997) results have suggested 
that clusters of well-fit genotypes that extend through- 
out genotype space are plausible. If this is so, biological 
populations are expected to evolve mainly within these 
clusters and consist most of the time of well-fit geno- 
types with fitnesses within some band. The metaphor 
of "holey" adaptive landscapes neglects the fitness differ- 
ences between genotypes in the cluster but these differ- 
ences are supposed to exist and should be apparent on a 
finer scale. If one applies a finer resolution, the movement 
along the cluster will be accompanied by slight increases 
or decreases in fitness. Evolution will proceed by fixa- 
tion of weakly selected alleles which can be advantageous, 
deleterious, over- and underdominant, or apparently neu- 
tral depending on the specific area of genotype space the 
population passes through. Smaller populations will pass 
faster through the areas of genotype space correspond- 
ing to fixation of slightly deleterious mutations whereas 
larger populations will pass faster through the areas cor- 
responding to fixation of (compensatory) slightly advan- 
tageous mutations. This pattern of molecular evolution, 
as predicted from the general properties of multidimen- 
sional adaptive landscapes, is similar to the patterns re- 
vealed by the methods of experimental molecular biology, 
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which form the empirical basis for the nearly neutral the- 
ory of molecular evolution (Ohta 1992). From general 
considerations, one should not expect complete symme- 
try of "real" adaptive landscapes which are supposed to 
have areas varying with respect to the "width" and con- 
centration of ridges of well-fit genotypes. Sexual popula- 
tions are expected to spend more time in areas of high 
concentration of well-fit genotypes (Peliti and Bastolla 
1994). One of the biological manifestations of this effect 
will be apparent reduction in the probability of harmful 
mutations, that is, evolution of genetic canalization (cf., 
Wagner 1996). The metaphor of holey adaptive land- 
scapes may be useful for thinking about these and other 
evolutionary problems. 
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Appendix 

Effects of mutation, migration and drift on the dynam- 
ics of the average genetic distances within and between 
subpopulations have been previously studied thoroughly 
(e.g., Watterson 1975; Li 1976; Slatkin 1987; Strobeck 
1987). What is left is to add reproductive isolation (that 
is selection) to the model. I will use the deterministic 
framework assuming that the population size N — > oo. 

The distribution of D w under rare- alleles and linkage 
equilibrium approximation. I will use the standard no- 
tations Ai and ai for alternative alleles at the i-locus 
(i = 1, . . . , L). Let pi be the frequency of allele A\ at the 
i-th locus, qi = 1—pi, and ip Wt i = 2piqi. Variable ip w ^ can 
be thought of as the probability that two randomly cho- 
sen individuals (sequences) from the same subpopulation 
are different at the i-th locus. Let d Wt i = (If — if) 2 be the 
genetic distance at the i-th locus between two randomly 
chosen individuals a and (3. Note that d w ^ = 1 with 
probability tp Wt i and d w j = with probability 1 — ip w ^. 
Because d w ^ is a Binomial random variable, its gener- 
ating function is "fd m ;(s) = tp w .i s + 1 — ipw,i, which can 
be approximated as e^™ '*( s_1 ) if tp Wti << 1 (rare-alleles 
approximation). Under approximate linkage equilibrium, 



the generating function of d w — ^2 dw,i is 

ld J s ) = n^.it-i) = ifr-'C'- 1 ) = e A.(-i) (Ai) 

where D w — Ylj ifrw.i- This shows that random variable 
d w has approximately Poisson distribution with parame- 
ter D w and, thus, 



P(d w = i) = exp(-D w )- 



(A2) 



Selection within an isolated population. Let w(d) be 
the expected number of fertile and viable offspring that 
can be produced by a pair of individuals different in d 
loci. The average fitness of the population is 



5>a)p(d=j). 



The dynamics of the general model of fertility selection 
and premating isolation in a haploid population consid- 
ered here are identical to that of a symmetric viability 
selection model for a diploid population with viabilities 
w(d) depending on the number of heterozygous loci d. 
Under approximate linkage equilibrium, changes in allele 
frequencies are described by Wright's equation 



Piqi d\nw 
~2 dp~' 



(A3) 



(Wright 1969). Using the equalities d\nw/dpi = 2(qi — 
Pi)d\nw/dipi and D w = Y^i^i) equation JA3|) can be 
rewritten as 

A s pi = spiqi(pi - qi), (A4a) 



with 



dhil 
~dT7 v 



(A4b) 



To describe the dynamics of allele frequencies one needs 
to know the mean fitness of the population. 

Truncation selection. This is a selection scheme ana- 
lyzed in detail in the main body of the paper. Here 



w{d) = 



1 for d < K, 
for d > K. 



(A5) 



Using the Poisson approximation (AI), the mean fitness 
is 



Wthreshold = eXp(-D w )- 
i=0 



T(K+l) ' 



where the last equality follows from equation (8.352) in 
Gradshteyn and Ryzhik (1994), resulting in s given by 
equation (pi]). To find equation (|l0|), one starts with 
(A4a) and proceeds using the fact that « 2(<^ — 
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Pi)Api and that D w =J2^i- 

Other selection schemes can be considered in a simi- 
lar way, and some of them result in relatively compact 
expressions for w and s. 

Linear selection. Here 



w(d) 



The mean fitness is 



1 - ad for d < K, 
for d > K. 



(A6) 



^linear — ^threshold &~ 



D W T(K,D U 
T(K) 



Quadratic selection. Here 
w(d) = 
The mean fitness is 

^ — ^linear 



I -ad- bd' 2 for d<K 
for d > K 



(A7) 



b[D w {D w -K) 
, D W (K + 1)T{K,D V 



-exp(-D w ) 



T(K) 

D%+ 2 H(2,K + 2,D W ) 



T(K + 2) 



where H is the hypergeometric function (Gradshteyn and 
Ryzhik 1994). 

Exponential selection. Here 



w(d) -- 
The mean fitness is 



exp(-ad) for d < K, 
for d > K. 



(A8) 



w = exp(—D w (l — e a )) 



T{K + l,D w e- a ) 
T{K + 1) 



I have not explored how assuming these selection schemes 
would affect the outcome of the dynamics. 

Stochastic transitions in an isolated population. 
Adding mutation results in equation 



Ap l = sp^iipi - qi) + (ify - p^, 



(A9) 



where [i is the rate of mutation (assumed to be equal for 
forward and backward mutations). Equation (|A9|) is sim- 
ilar to the classical equation describing underdominant 
selection on a single locus in a diploid population. This 
allows one to use Lande's results (1979; see also Hedrick 
1981; Walsh 1982; Barton and Rouhani 1987) to find the 
rate of stochastic divergence. This rate is twice the ex- 
pected number, vN, of new mutations in a population 
times the probability that a given one will be fixed, U. 
Using the diffusion approximation, U is defined by equa- 



tions (la) and (2) in (Lande 1979). Lande used some 
approximations to evaluate U. However, the integrals in 
his equation (la) can be found exactly resulting in 



U 



erf 


VS(1- 


N > 


erf 


Vs 





(A10) 



where S = Ns/2 (Walsh 1983). Expanding in a Taylor 
s erie s under the assumption that 1/N << 1 results in 
(13b) which is equivalent to Lande's (1979) formula. The 
difference between Lande's approximate formula and the 
exact equation ( A10 ) is negligible. 

The distribution of Db under rare-alleles and linkage 
equilibrium approximation. Let us consider two subpop- 
ulations. Let pi and Pi be the frequencies of allele Ai 
in the first and second subpopulations, respectively. The 
genetic distance db,i at the i-th locus between two ran- 
domly chosen sequences from two different subpopula- 
tions is a Binomial variable taking values 1 and with 
probabilities ipb.i = PiQi + qiPi and 1 — ipb,i, respectively 
(qi = 1—pi, Qi = I — Pi). I will assume that genetic varia- 
tion within each subpopulation is low so that ipb,i is close 
to either or 1. Let Si — db,i if ipb,i ~ and Si = 1 — db.i if 
tpb,i ~ 1- The genetic distance between individuals a and 
P can be represented as db — k — J^i $i + TI2 ^» wnere the 
first sum is over k loci at which ipb,i ~ 1 an d the second 
sum is over L — k loci at which ipb,i ~ 0. Using the as- 
sumption of linkage equilibrium, the generating function 
of db becomes 



7* 



(s) =^{ s fc -E 1 ^+E 2 ^} 



= s k Ii k ; =1 e-^ {s - 1) IV^ =k+1 e^ {s - 1) 



s _ e D b (s-l) 



gfc(s-l) 

where <f>i is the expectation of Si and Db is the expectation 
of db. Using the properties of generating functions, the 
distribution of db is 







if i < k, 



(All) 



With fitness function ( |A5|) , the probability that two ran- 
domly chosen individuals from different subpopulations 
are not reproductively isolated is 



K 



b = t) 



i=0 



T{K - k + 1, D b - k) 
T(K-k + l) ' 



(A12) 



if k < G and Wb = if k > K . 

Deterministic dynamics in a subdivided population. 
With no reproductive isolation and with equal forward 
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and backward migration rates and equal population sizes, 
the change in pi due to migration is A m pi = m(Pi ~ Pi). 
The corresponding change in D w is A m D w — 2m(Db — 
D w ). With reproductive isolation and given that Df, > 
D w , individuals migrating from other subpopulations will 
have reduced probability of mating. Let W w and w b be 
the expected numbers of fertile and viable offspring that 
can be produced as a result of within and between sub- 
populations encounters. For simplicity I will omit the in- 
dex specifying the locus under consideration. With equal 
population sizes and migration rates, the change in the 
allele frequency due to migration becomes 



In the former k loci, the dynamics of ip w are described as 



A m p = m e (P - p) , 
where the "effective" migration rate is 



(Al3a) 



(Al3b) 



(compare with the models with migration between popu- 
lations of unequal size where the effective migration rate 
is m times the ratio of the population sizes, e.g. Gavrilets 
1996). The corresponding change in D w is 



A m D w 



2m e (D b -D w ). 



(A13c) 



Changes A m pi can be thought of as changes in allele fre- 
quencies brought about by selection between groups of in- 
dividuals (migrants and residents) whe reas the first term 
in the right-hand side of equation (A9) can be thought of 
as the change in p brought about by individual selection. 



The dynamics of allele frequencies at a specific locus 
under the joint action of selection, mutation and migra- 
tion arc described by 

Ap = spq(p - q) + m e (P — p) + fJ,(q - p), (Al4a) 
AP = ~sPQ(P -Q) + m e (p-P)+ n(Q - PjAUb) 

With s = $ = c onst and m e = const and m e < s/6, 



dynamic system (AAA) has two types of stable equilibria: 
mutation-selection balance equilibria with p w P, pq w 
/V s ! i>w ~ w 2/i/s and migration-selection equi- 

libria with p s» Q, pq ss u/s + m/s, tp w ps 2/i/s + 
2m/ s, ipj, ps 1 — 2/x/s — 2m/ s. I assume that in the deter- 
ministic limit, k out of L loci evolve towards migration- 
selection balance equilibria whereas the remaining L — k 
loci evolve towards mutation-selection balance equilibria. 
In the latter L — k loci, the dynamics of ip w and ip b are 
approximated by equations 

Aip w = -sip w + 2(i + 2m e (ip b - ipw), (A15a) 
Aip b = -sipb + 2/i + 2m e (i/v, - ip b ). (Al5b) 



before by ( A15a) whereas the dynamics of ipb are approx- 
imated by equation 

Aip h = s{l-ip h )-2ix + 2m e {ip w - ip h ). (A16) 



Selection always reduces ip 



increases it (see equation Al5a ). 
tion have the same effects on ipb 
towa rd mu tation-selection balance equilibria (see equa- 
However, for the loci evolving toward 



whereas mutation always 
Selection and muta- 
for the loci evolving 



A15b ) 



tion 

migration-selection balance equilibria, selection inc rease s 
ipb whereas mutation decreases it (see equation A16 ). 
Summing up over all loci, one finds equations ([l4|) of the 
main text. Equations (^) are derived in a similar way as- 
suming that the allele frequencies in the main population 
do not change. 

Stochastic transitions in a subdivided population. In 
a subdivided population, migration tends to reduce ge- 
netic differentiation. Given that migration is sufficiently 
strong relative to selection, the same allele will be close 
to fixation in both subpopulations. If, by a chance, an 
alternative allele approaches fixation in one of the sub- 
populations creating significant differentiation at a given 
locus, such differentiation will be quickly lost. The num- 
ber of loci k at which alternative alleles are close to fix- 
ation in different subpopulations will be close to zero on 
average. However, if migration is relatively weak, then 
the differentiation created by random genetic drift will 
not be lost quickly and actually can even accumulate. 
Let us consider a locus at which initially the same allele 
is close to fixation in both subpopulations (that is both 
p ps and P rj 0). Neglecting the changes in P, the de- 
terministic change in p due to selection and migration is 
approximately 



Ap = spq(p — q) — m e p. 



(A17) 



Lande (1979; see also Barton and Rouhani 1987) has 
shown that the rate at which allele A becomes close to 
fixation in the first subpopulation while its frequency is 
about zero in the second population is approximately 
2~Nm a tirngg the rate of fixation in the absence of im- 
migration. Assuming that alleles A are brought about 
by mutation at rate \i and summing up over L — k loci, 
one finds the first term in the right-hand side of equation 
(P~r|) . Once alternative alleles are close to fixation in dif- 
ferent subpopulations, random drift can remove genetic 
differentiation. Let us consider a locus at which initially 
p w but Pal. Neglecting the changes in P, the deter- 
ministic change in p due to selection and immigration is 
approximately 



Ap = spq(p - q) - m e (1 - p). 



(A18) 



February 5, 2008 



Page 16 



Using Barton and Rouhani's (1987) method one finds that 
the rate at which allele A becomes close to fixation in 
both subpopulations is approximately (e/2) Nmc times the 
rate of fixation in the absence of immigration. Assuming 
that alleles A are brought about by migration at rate m e 
and summing up over k loci, one finds the second term in 
the right-hand side of equation (p7|). 

Stochastic divergence with local adaptation. Let us as- 
sume that the allele under consideration is favorable in 
a given environment with selective advantage sla- The 
change in this allele frequency as defined by the joint ac- 
tion of selection induced by reproductive isolation and 
selection for local adaptation is 

A s p = spq(p - q) + s LA pq. (A19) 

This equation is identical to the one describing meiotic 
drive in the Appendix of Walsh (1982). Following Walsh, 
the fixation probability is 



erf 


VS(l-a) 


- erf 


Vs(i -<*-#)" 


erf 


y/S{l-a) 


+ erf 


VS{l + a) 





(A20) 



where S — Ns/2 and a = Sla/s. Expanding the nu- 
merator in a Taylor series under the assumption that 
1/N << 1 and multiplying the results by the expected 
number of mutants, vN, results in the relative rate of 
fixation given by equation (|l9|). 

Genetic distance ( pZ| j and some other models. Genetic 
distance (1) is recovered by assuming that G is an iden- 
tity matrix. Assuming that G is a diagonal matrix with 
non-equal diagonal elements is a simple way to introduce 
non-equivalence of loci. The case when the probability of 
mating depends on the difference in a quantitative trait 
can be treated within the same framework. Let c; be the 
contribution of the i-th locus to a quantitative trait z. 
Neglecting microenvironmental effects, the trait values for 
individuals a and (3 are x a — c%lf and z' 3 = Cjif , re- 
spectively. The square of the difference of z a and z' 3 is re- 
covered from equation ( ^p|) by assuming that G^j — c^Cj. 
A common way to model sexual selection is to assume 
that the probability of mating between a male and a fe- 
male depends on the difference in a female phenotypic 
trait, Zf, and a male phenotypic trait, z m , which are 
controlled by two different sets of loci (e.g., Lande 1981; 
Kirkpatrick 1982; Nei at al. 1983; Wu 1985; Turner and 
Burrows 1995). Let z m = c^h and zt = c{k where 
the sums are taken over the corresponding sets of loci. 
The value (z m — zj) 2 is recovered from (po| ) by assuming 
that matrix G has a block form 

_ / Q G s \ 
\ G s J ' 



The diagonal L m x L m and Lf x Lf zero matrices cor- 
respond to the interactions within the set of L m genes 
controlling the male trait and within the set of Lf genes 
controlling the female trait (L = Lf+L m ), and matrix G s 
describing the interactions between the two sets of genes 
has elements Gjj = c{c 1 J l . 
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Figure 6: Dynamics of speciation in a subdivided popula- 
tion. Unless specified otherwise, K = 20, v = 0.0384, n = 2. 
a. Effects of migration rate. Stronger migration, m = 
0.01 (dashed lines; no speciation), and weaker migration, 
m = 0.001 (solid lines; speciation). Total population size 
Nt = 200. b. Effects of mutation rate. Weaker mutation, 
v — 0.0384 (dashed lines; no speciation), and stronger muta- 
tion, v — 5 x 0.0768 (solid lines; speciation). Other param- 
eters: Nt = 400, m = 0.005. c. Effect of population sub- 
division, n — 2 subpopulations (dashed lines; no speciation) 
and n = 4 subpopulations (solid lines; speciation). Other pa- 
rameters: Nt = 800, m = 0.0033. Dashed lines represent Df, 
(top line), D w (middle line) and k (bottom line), respectively. 
During the first 1000 generations there are no restrictions on 
migration. 
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Figure 7: Dynamics of speciation in a peripheral population, 
a. Speciation with m = 0.001; N = 100; K = 20; v = 0.0384 
(cf. Fig. 7a). b. Speciation with m = 0.005; N = 200; if = 
20; i; = 5 x 0.0768 (cf. Fig. 7b). 
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Figure 8: Relative rate of fixation in the case with local 
adaptation. 



